Méthode de Simpson

La méthode des rectangles approche chaque portion de courbe par une fonction constante, c’est-à-dire une fonction polynômiale de degré au plus 0. La méthode des trapèzes, quant à elle, approche chaque portion de courbe par une fonction affine, éventuellement constante. En d’autres termes, on approche la fonction à l’aide d’une fonction polynomiale de degré au plus 1. On peut alors naturellement se demander si ces méthodes peuvent s’étendre à des fonctions de degré supérieur… Et en effet, elles le peuvent !

La méthode de Simpson approche chaque portion de courbe par un arc de parabole ou un segment de droite, selon les cas. La fonction étudiée est donc approchée par une fonction polynômiale de degré au plus 2.

Exercice
On considère un intervalle  \([s;t]\) et on note  \(m=\dfrac{s+t}{2}\) .
On considère la fonction  \(P:x\mapsto f(s)\dfrac{(x-t)(x-m)}{(s-t)(s-m)} +  f(t)\dfrac{(x-s)(x-m)}{(t-s)(t-m)} +  f(m)\dfrac{(x-t)(x-s)}{(m-t)(m-s)}\)

1. La fonction  \(P\) est-elle polynômiale ? Quel est son degré ?

2. Calculer  \(P(t)\) ,  \(P(s)\)  et  \(P(m)\) .

On utilise alors la fonction  \(P\) définie ci-dessus pour approcher la fonction  \(f\) sur l'intervalle \([s;t]\) .
L'intégrale de \(f\)  sur \([s;t]\)  peut alors être approchée par l'intégrale de `P` sur ce même intervalle.
On peut par ailleurs montrer que \(\displaystyle\int_s^t P(x)dx = \dfrac{t-s}{6} \times (f(s) + 4f(m) + f(t))\) .
Là encore, on peut l’implémenter en Python pour se rendre compte de la vitesse de convergence de cette méthode, qui surpasse les vitesses de convergence des deux méthodes précédentes. Dans le cas suivant, on utilise simplement 3 subdivisions de l'intervalle \([1;4]\) .

from math import log, sqrt

def fonction(x) :
    return  2*x**2 *log(x) + x - 2*sqrt(x)

def simpson(f, a, b, n) :
    '''
    Approche l'intégrale de f entre a et b en subdivisant l'intervalle [a,b]
    en n sous-intervalles et en utilisant la méthode de Simpson
    ----------
    Entrées :
    f : une fonction continue dont on souhaite calculer l'intégrale
    a : un réel, borne gauche de l'intervalle d'intégration
    b : un réel, borne droite de l'intervalle d'intégration
    n : un entier, nombre de subdivisions de l'intervalle [a,b]
    ----------
    Sorties :
    I : un réel, valeur approchée de l'intégrale de f
    '''
    I = 0
    dx = (b-a) / n
    x = a
    for i in range(n) :
        I += (f(x) + f(x + dx) + 4 * f(x + dx/2)) * dx / 6
        x += dx
    return I

print(simpson(fonction, 1, 4, 3))
#Valeur réelle de l'intégrale : 128ln(64)-95/6, soit environ 43.315